#!/usr/bin/perl
# 2011/08/04 Overlapしている数が有意かどうかを検定する

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);
use List::Util qw(shuffle);
$| = 0;


if((@ARGV != 6 and @ARGV != 5) or $ARGV[0] eq '--help'){
	die "Usage : $0 -a [bed file #1] -b [bed file #2] -o [output image name]\n";
}

my $search_width = 300;

my %opt;
getopts("a:b:o:", \%opt);
my $FILE_A = $opt{a};
my $FILE_B = $opt{b};
my $FILE_image = $opt{o};

#---------------------------------------
# bed fileの読み込み
#---------------------------------------
my %all;
my %sample = &readBed($FILE_A, \%all);
my %target = &readBed($FILE_B, \%all);
sub readBed{
	my ($file, $ref) = @_;
	my %hash;
	my $number = 0;
	my $fh = IO::File->new($file) or die "cannot open $file: $!";
	while($_ = $fh->getline()){
		s/\r?\n//;
		if(m/^track/){
			next;
		}

		$number++;
		my ($chr, $start, $end, @other) = split /\t/;
		my $start_c = int($start / 1000) * 1000 + 500;
		my $end_c = int($end / 1000) * 1000 + 500;

		for(my $i = $start_c; $i <= $end_c; $i += 1000){
			$hash{$chr}{$i}++;
			$ref->{$chr}{$i} = 1;
		}
	}
	$fh->close();
	return %hash;
}

#---------------------------------------
# ファイルを出力
#---------------------------------------
my %MAX = (
	'chr1'	=>	5672121,
	'chr2'	=>	4615470,
	'chr3'	=>	2493767,
);

my $tmp = 'tmp.txt';
my $fh_out = IO::File->new($tmp, 'w') or die "cannot write $tmp: $!";
$fh_out->print("num\tsam\ttarget\n");

my $number = 0;
foreach my $chr(qw(chr1 chr2 chr3)){
	for(my $c = 500; $c <= int($MAX{$chr} / 1000) * 1000 + 500; $c += 1000){
		$number++;
		$fh_out->printf("%d\t%d\t%d\n", $number, exists $sample{$chr}{$c} ? $sample{$chr}{$c} : 0, exists $target{$chr}{$c} ? $target{$chr}{$c} : 0);
	}
}
$fh_out->close();

my $Rcommand_p = <<END;
R --vanilla --slave <<EOF
	data <- read.table("$tmp", header=T, sep="\t")
	attach(data)

	n = 10000

	sample_num = length(sam[sam > 0])
	common_num = length(sam[sam > 0 & target > 0])

	empty = rep(1, length(sam))

	count = 0
	distribution = c()
	for(i in 1:n){
		random = sample(num, sample_num, replace=FALSE)
		com = sum(target[random] > 0)
		distribution = c(distribution, com)
		if(com > common_num){
			count = count + 1
		}
	}

	xmax = floor(common_num * 1.1)
	if(max(distribution) > xmax){
		xmax = max(distribution)
	}

	pvalue = count / n
#	jpeg(file="$FILE_image", width=700, height=400)
	pdf(file="$FILE_image")
	hist(distribution, xlim=c(0,xmax), main="common peak number distribution")
	abline(v=common_num, col="red")
	legend("topright", lty=1, lwd=2, col="red", legend=paste("pvalue = ", sprintf("%.3f", pvalue), sep=""))
#	graphics.off()
	dev.off()
	cat(pvalue)

	#write.table(distribution, file="out.txt", quote=F, sep="\t", row.names=F, col.names=F)
EOF
END

my $pvalue = `$Rcommand_p`;
print "p-value = $pvalue\n";
unlink $tmp;











